rm(list = ls())
#set working directory
setwd("Replication Data")
#setwd("~/Dropbox/War and Hierarchies/ISQ_Final_Submission/Replication Data")

## Load the R packages
library(readstata13)
library(arm)
library(ggplot2)
library(ggthemes)
library(viridis)
library(ggridges)
library(purrr)
library(dplyr)
library(readxl)
library(readr)
library(arm)
library(plm)
library(haven)
library(gridExtra)
library(grid)
library(texreg)
#load country-level data for analysis
load("data_cty.RData")

##### load R function to run the model
source("RFunctions.R")

############### fit fixed-effect model with robust standar errors by group
# return as a ggplot objective

###########################################################################
## Fig 1: war and changes in civil war libterties
###########################################################################

# Fig 1c: women's civil liberties
#############################################
#interstate: women's civil liberties 
inter_war_p <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                          'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                          'f10s11_civilliberty'),
                ivs = c('inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                        "s_lpec", "l_lpec", "yrs"),
                modname = c('womenCL0', 'womenCL1', 'womenCL2', 'womenCL3', 'womenCL4',
                            'womenCL5', 'womenCL10'),
                n_sim = 1000,
                var = "inter_warDummy",
                data = data_cty)
#add labels
p1 <- inter_war_p$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year",
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Interstate War",
            x = "Average marginal effects")

#check result with Stata output: confirmed 
inter_war_p$Results[[1]]
summary(inter_war_p$Fits[[1]])
#get individual intercept
fixef(inter_war_p$Fits[[1]])
mean(fixef(inter_war_p$Fits[[1]]))


#### for Intrastate war: women's civil liberties 
intra_war_p <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                  'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                  'f10s11_civilliberty'),
                          ivs = c('inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                                  "s_lpec", "l_lpec",'yrs'),
                          modname = c('womenCL0', 'womenCL1', 'womenCL2', 'womenCL3', 'womenCL4',
                                      'womenCL5', 'womenCL10'),
                          n_sim = 1000,
                          var = "intra_warDummy",
                          data = data_cty)
## add labels
p2 <- intra_war_p$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
     labs(title = "Intrastate War",
       x = "Average marginal effects")
## save and arrange two figures together

f <-grid.arrange(p1, p2, nrow = 1, 
             top = textGrob("War and Changes in Gender Civil-liberty Equality",
                            gp = gpar(face = "bold", fontface = 14, fontsize = 14),
                            hjust = -0.2,x = 0))
#dev.off()
ggsave("figures/Figure_1_bottom.eps", plot = f, width = 6.5, height = 3.5, dpi = 400, units = "in")


#### appendix fig 1c: make coefficient plot for women's civil liberties (Appendix Fig 1)
#shape palette can deal with a maximum of 6 
# The order of the labels will be the order of the input
p <- coef_rb_plot(inter_war_p$Results, vars = c('inter_warDummy',
                                                'intra_warDummy',
                                                "s_v2x_polyarchy", 
                                                "l_v2x_polyarchy",
                                                "s_lpec", 
                                                "l_lpec",
                                                "yrs")) 
#add lables 
p <- p + scale_x_discrete(labels=c("Interstate War", "Intrastate war",
                              "Polyarchy score (change)", "Polyarchy score (lag)",
                              "Energy (change)", "Energy(lag)", "Year")) + 
       ggtitle("Coefficient Plot: War and Changes in Gender Civil-liberty Equality") + 
       scale_color_discrete(name = "", 
                            labels = c("Current", "1-year",  "2-year",
                                          "3-year", "4-year", "5-year",
                                    "10-year"))
ggsave("figures/Appendix_Fig_1c.pdf", plot = p, width = 8, height = 6)
#############################################################################



# Fig 1a: war and ethnic civil liberties
# interstate: ethnic civil liberties 

inter_war_q <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                  'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                  'f10s11_v2clsocgrp'),
                          ivs = c('inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                                  "s_lpec", "l_lpec", "yrs"),
                          modname = c('ethnicityCL0', 'ethnicityCL1', 'ethnicityCL2', 'ethnicityCL3', 'ethnicityCL4',
                                      'ethnicityCL5', 'ethnicityCL10'),
                          n_sim = 1000,
                          var = "inter_warDummy",
                          data = data_cty)
# add labels
p3 <- inter_war_q$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
                       labs(title = "Interstate War",
                            x = "Average marginal effects")

# intrastate: ethnic civil liberties 

intra_war_q <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                  'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                  'f10s11_v2clsocgrp'),
                          ivs = c('inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                                  "s_lpec", "l_lpec", "yrs"),
                          modname = c('ethnicityCL0', 'ethnicityCL1', 'ethnicityCL2', 'ethnicityCL3', 'ethnicityCL4',
                                      'ethnicityCL5', 'ethnicityCL10'),
                          n_sim = 1000,
                          var = "intra_warDummy",
                          data = data_cty)

p4 <- intra_war_q$plot + scale_y_discrete(labels = c("10-year", "5-year",
                                                     "4-year", "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Intrastate War",
       x = "Average marginal effects")

## save and arrange two figures together
#pdf(file = "figures/RR/marg_FE_ethnicCL.pdf", width = 8, height = 6)
f <- grid.arrange(p3, p4, nrow = 1, 
             top = textGrob("War and Changes in Ethnic Civil-liberty Equality",
                            gp = gpar(face = "bold", fontface = 14, fontsize = 14),
                            hjust = -0.2,x = 0))
#dev.off()
ggsave("figures/Figure_1_top.eps", plot = f, width = 6.5, height = 3.5, dpi = 400, units = "in")



## appendix fig 1a
#make coefficient plot
q <- coef_rb_plot(inter_war_q$Results, vars = c('inter_warDummy',
                                                'intra_warDummy',
                                                "s_v2x_polyarchy",
                                                "l_v2x_polyarchy",
                                                "s_lpec", 
                                                "l_lpec",
                                                "yrs")) 

q <- q + scale_x_discrete(labels=c("Interstate War", "Intrastate war",
                              "Polyarchy score (change)", "Polyarchy score (lag)",
                              "Energy (change)", "Energy(lag)", "Year")) + 
         ggtitle("Coefficient Plot: War and Changes in Ethnic Civil-liberty Equality") + 
         scale_color_discrete(name = "", 
                              labels = c("Current", "1-year",  "2-year",
                                         "3-year", "4-year", "5-year",
                                         "10-year"))
ggsave("figures/Appendix_Fig_1a.pdf", plot = q, width = 8, height = 6)

#################################################################
############################################################################################
# Fig 1 b: war and ethnic civil liberties when women's empowerment increase

dvs <- c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
        'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
        'f10s11_v2clsocgrp')
ivs <- c('inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy",
         "l_v2x_polyarchy", 
        "s_lpec", "l_lpec", "yrs")
modname <- c('ethnicityCL0_int', 'ethnicityCL1_int', 'ethnicityCL2_int', 'ethnicityCL3_int', 'ethnicityCL4_int',
            'ethnicityCL5_int', 'ethnicityCL10_int')
conditions <- c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                'f10s11_civilliberty')

n_sim = 1000
var <- c('inter_warDummy' , 'intra_warDummy')

#empty list to store marginal results       
marginal <- list()     
#empty list to store coefficents
Results <- list()
# empty list to store model fit
Fits <- list()
fig <- list()

# for loop to run all models

for (j in seq_along(var)) {

for (i in seq_along(modname)){
       f <- as.formula(paste(dvs[i], paste(ivs, collapse=" + "), sep=" ~ "))
        
       data <- data_cty %>% 
           dplyr::filter(get(conditions[i]) > 0 & !is.na(get(conditions[i])))
       # fit a fixed effect model at country-year-level
       fit <- plm(f, data = data, model = "within", index = c("ccode", "year"))
       rb_coef <- coeftest(fit, vcov=vcovHC(fit, type="sss", cluster="group"))
       # store results
       Results[[modname[i]]] <- rb_coef
       #extract the variance-covariance matrix
       vc <- vcovHC(fit, type="sss", cluster="group")
       #extract the coefficient 
       pars <- rb_coef[, "Estimate"]
       
       #extrat the model data
       mydata2 <- model.matrix(f, data)
       #remove constant 1
       mydata2 <- as.data.frame(mydata2[,c(2: (length(pars) +1))])# remove constant 

       ### simulation based on observed values
       # simulations based on coefficient and variance-covariance
       set.seed(321)
       draw <- mvrnorm(n_sim, data.matrix(pars), data.matrix(vc))
       #copy dataset for marginal effect
       data0 <- data1 <- mydata2
       #set the scenario for marginal effect
       data0[, var[j]] <- 0
       data1[, var[j]] <- 1
       
       #empty containers
       X1 <- as.matrix(data0)
       X2 <- as.matrix(data1)
       #store marginal effect in the list
       marginal[[modname[i]]] <- apply(apply(X2, 1, function (x) draw %*% x) - 
                                              apply(X1, 1, function (x) draw %*% x), 1, mean) #1 indicates row, 2= columns
       
       
       }

#convert the list to a data frame for plotting
df_plot <- map(marginal, data.frame) %>%
       map2_df(., names(.), ~mutate(.x, type = .y))
#rename first variable
names(df_plot)[1] <- "value"

##change the order of the variable
df_plot$type <- factor(df_plot$type, levels = rev(modname))

#### Plotting
fig[[var[j]]] <- ggplot(df_plot, aes(x = value, y = type, height=..density.., fill = type)) +
       geom_density_ridges(col = "grey70", scale = 1.2, show.legend = F) +
       scale_fill_viridis(discrete = TRUE) +
       geom_vline(xintercept = 0, colour = gray(1/2), lty = 2) +
       theme_ridges(font_size = 13, grid = F, center_axis_labels = T) +
       theme(axis.title.y = element_blank(), 
             text = element_text(size=13),
             plot.title = element_text(hjust = .5, size = 15, face = "bold")) 
}

# make figures
# add labels
p1 <- fig$inter_warDummy + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Interstate War",
            x = "Average marginal effects")

p2 <- fig$intra_warDummy + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                       "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Intrastate War",
            x = "Average marginal effects")
## save and arrange two figures together
#pdf(file = "figures/RR/marg_FE_ethnicCL_womenCL.pdf", width = 8, height = 6)
f <- grid.arrange(p1, p2, nrow = 1, 
             top = textGrob("War and Changes in Ethnic Civil-liberty Equality\n When Gender Civil-liberty Equality Improves",
                            gp = gpar(face = "bold", fontface = 14, fontsize = 14),
                            hjust = -.4,x = 0))
#dev.off()
ggsave("figures/Figure_1_middle.eps", plot = f, width = 6.5, height = 3.5, dpi = 400, units = "in")


# make a coefficient plot: Appendix fig 1b
q2 <- coef_rb_plot(Results, vars = c('inter_warDummy',
                                      'intra_warDummy',
                                     "s_v2x_polyarchy",
                                     "l_v2x_polyarchy",
                                       "s_lpec", 
                                         "l_lpec",
                                          "yrs")) 

q2 <- q2 + scale_x_discrete(labels=c("Interstate War", "Intrastate war",
                                   "Polyarchy score (change)", "Polyarchy score (lag)",
                                   "Energy (change)", "Energy(lag)", "Year")) + 
       ggtitle("Coefficient Plot: War and Changes in Ethnic Civil-liberty Equality\n When Gender Civil-liberty Equality Improves") + 
       scale_color_discrete(name = "", 
                            labels = c("Current", "1-year",  "2-year",
                                       "3-year", "4-year", "5-year",
                                       "10-year"))
ggsave("figures/Appendix_Fig_1b.pdf", plot = q2, width = 10, height = 7)




####################################################################################
######################### Fig 2
## Figure 2a: interstate war outcomes for group civlib 
recent_loss_inter <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                        'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                        'f10s11_v2clsocgrp'),
                                ivs = c('new_inter' , 'ongoing_inter',"recent_win_inter","recent_loss_inter","recent_draw_inter",
                                        "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                                modname = c('ethnicCL_inter0', 'ethnicCL_inter1', 'ethnicCL_inter2', 'ethnicCL_inter3', 'ethnicCL_inter4',
                                            'ethnicCL_inter5', 'ethnicCL_inter10'),
                                n_sim = 1000,
                                var = "recent_loss_inter",
                                data = data_cty)
# add labels
p2 <- recent_loss_inter$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                           "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Recent Interstate War Losses and Changes in Ethnic Civil-liberty Equality",
       x = "Average marginal effects")
ggsave("figures/Figure_2a.eps", plot = p2, width = 6.5, height = 3.5, dpi = 400, units = "in")


## coefficient plot: Figure A9-a

q5 <- coef_rb_plot(recent_loss_inter$Results, vars = c('new_inter' , 'ongoing_inter',"recent_win_inter","recent_loss_inter","recent_draw_inter",
                                                       "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs")) 

q5 <- q5 + scale_x_discrete(labels=c("New interstate war", "Ongoing interstate war", "Recent interstate war win", "Recent interstate war loss",
                                     "Recent interstate war draw","Polyarchy score (change)", "Polyarchy score (lag)",
                                     "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: Interstate War Outcomes and Changes in Ethnic Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_9a.pdf", plot = q5, width = 11, height = 9)


### Figure 2b: intrastate war outcome for group 
recent_loss_intra <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                        'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                        'f10s11_v2clsocgrp'),
                                ivs = c('new_intra' , 'ongoing_intra',"recent_win_intra","recent_loss_intra","recent_draw_intra",
                                        "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                                modname = c('ethnicCL_intra0', 'ethnicCL_intra1', 'ethnicCL_intra2', 'ethnicCL_intra3', 'ethnicCL_intra4',
                                            'ethnicCL_intra5', 'ethnicCL_intra10'),
                                n_sim = 1000,
                                var = "recent_loss_intra",
                                data = data_cty)
# add labels
p2 <- recent_loss_intra$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                           "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Recent Intrastate War Losses and Changes in Ethnic Civil-liberty Equality",
       x = "Average marginal effects")
ggsave("figures/Figure_2b.eps", plot = p2, width = 6.5, height = 3.5, dpi = 400, units = "in")


## ## coefficient plot: Figure A9-b
q6 <- coef_rb_plot(recent_loss_intra$Results, vars = c('new_intra' , 'ongoing_intra',"recent_win_intra","recent_loss_intra","recent_draw_intra",
                                                       "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs")) 

q6 <- q6 + scale_x_discrete(labels=c("New intrastate war", "Ongoing intrastate war", "Recent intrastate war win", "Recent intrastate war loss",
                                     "Recent intrastate war draw","Polyarchy score (change)", "Polyarchy score (lag)",
                                     "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: Intrastate War Outcomes and Changes in Ethnic Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_9b.pdf", plot = q6, width = 11, height = 9)
############################################################################################

## Figure 2c: interstate war loss and women's civil 
recent_loss_inter <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                        'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                        'f10s11_civilliberty'),
                                ivs = c('new_inter' , 'ongoing_inter',"recent_win_inter","recent_loss_inter","recent_draw_inter",
                                        "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                                modname = c('womenCL_inter0', 'womenCL_inter1', 'womenCL_inter2', 'womenCL_inter3', 'womenCL_inter4',
                                            'womenCL_inter5', 'womenCL_inter10'),
                                n_sim = 1000,
                                var = "recent_loss_inter",
                                data = data_cty)
# add labels
p2 <- recent_loss_inter$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                           "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Recent Interstate War Losses and Changes in Gender Civil-liberty Equality",
       x = "Average marginal effects")
ggsave("figures/Figure_2c.eps", plot = p2, width = 6.5, height = 3.5, dpi = 400, units = "in")

## coefficient plot: Figure A9-c

q7 <- coef_rb_plot(recent_loss_inter$Results, vars = c('new_inter' , 'ongoing_inter',"recent_win_inter","recent_loss_inter","recent_draw_inter",
                                                       "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs")) 

q7 <- q7 + scale_x_discrete(labels=c("New interstate war", "Ongoing interstate war", "Recent interstate war win", "Recent interstate war loss",
                                     "Recent interstate war draw","Polyarchy score (change)", "Polyarchy score (lag)",
                                     "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: Interstate War Outcomes and Changes in Gender Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_9c.pdf", plot = q7, width = 11, height = 9)

### Figure 2d: intrastate war loss and women's civil
recent_loss_intra <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                        'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                        'f10s11_civilliberty'),
                                ivs = c('new_intra' , 'ongoing_intra',"recent_win_intra","recent_loss_intra","recent_draw_intra",
                                        "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                                modname = c('womenCL_intra0', 'womenCL_intra1', 'womenCL_intra2', 'womenCL_intra3', 'womenCL_intra4',
                                            'womenCL_intra5', 'womenCL_intra10'),
                                n_sim = 1000,
                                var = "recent_loss_intra",
                                data = data_cty)
# add labels
p2 <- recent_loss_intra$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                           "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Recent Intrastate War Losses and Changes in Gender Civil-liberty Equality",
       x = "Average marginal effects")
ggsave("figures/Figure_2d.eps", plot = p2, width = 6.5, height = 3.5, dpi = 400, units = "in")

##  coefficient plot: Figure A9-d

q8 <- coef_rb_plot(recent_loss_intra$Results, vars = c('new_intra' , 'ongoing_intra',"recent_win_intra","recent_loss_intra","recent_draw_intra",
                                                       "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs")) 

q8 <- q8 + scale_x_discrete(labels=c("New intrastate war", "Ongoing intrastate war", "Recent intrastate war win", "Recent intrastate war loss",
                                     "Recent intrastate war draw","Polyarchy score (change)", "Polyarchy score (lag)",
                                     "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: Intrastate War Outcomes and Changes in Gender Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_9d.pdf", plot = q8, width = 11, height = 9)



###################################################
## Figure 3 war, regime change and change in ethnic civil war liberties
###########################################################################################################
###########################################################################################################
#regime change
regime <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                             'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                             'f10s11_v2clsocgrp'),
                     ivs = c('irregular_dummy', 'inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                             "s_lpec", "l_lpec", "yrs"),
                     modname = c('ethnicCL0', 'ethnicCL1', 'ethnicCL2', 'ethnicCL3', 'ethnicCL4',
                                 'ethnicCL5', 'ethnicCL10'),
                     n_sim = 1000,
                     var = "irregular_dummy",
                     data = data_cty)
# add labels
p_regime <- regime$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                      "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Regime Change",
       x = "Average marginal effects")

# intrastate
intra_war <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                'f10s11_v2clsocgrp'),
                        ivs = c('irregular_dummy', 'inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                                "s_lpec", "l_lpec", "yrs"),
                        modname = c('ethnicCL0', 'ethnicCL1', 'ethnicCL2', 'ethnicCL3', 'ethnicCL4',
                                    'ethnicCL5', 'ethnicCL10'),
                        n_sim = 1000,
                        var = "intra_warDummy",
                        data = data_cty)

p_intra_war <- intra_war$plot + scale_y_discrete(labels = c("10-year", "5-year",
                                                            "4-year", "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Intrastate War",
       x = "Average marginal effects")

#Interstate
inter_war <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                'f10s11_v2clsocgrp'),
                        ivs = c('irregular_dummy', 'inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                                "s_lpec", "l_lpec", "yrs"),
                        modname = c('ethnicCL0', 'ethnicCL1', 'ethnicCL2', 'ethnicCL3', 'ethnicCL4',
                                    'ethnicCL5', 'ethnicCL10'),
                        n_sim = 1000,
                        var = "inter_warDummy",
                        data = data_cty)

p_inter_war <- inter_war$plot + scale_y_discrete(labels = c("10-year", "5-year",
                                                            "4-year", "3-year", "2-year", "1-year", "Current")) + 
  labs(title = "Interstate War",
       x = "Average marginal effects")


## save and arrange two figures together
#pdf(file = "figures/RR/marg_FE_ethnicCL_regime.pdf", width = 9, height = 6)
f <- grid.arrange(p_regime, p_inter_war, p_intra_war, nrow = 1, 
             top = textGrob("War, Regime Change, and Changes in Ethnic Civil-liberty Equality",
                            gp = gpar(face = "bold", fontface = 12, fontsize = 12),
                            hjust = -0.2,x = 0))
ggsave("figures/Figure_3.eps", plot = f, width = 6.5, height = 3.5, dpi = 400, units = "in")
##dev.off()

#Appendix Figure A6
q10 <- coef_rb_plot(inter_war$Results, vars = c('irregular_dummy', 'inter_warDummy' , 'intra_warDummy',"s_v2x_polyarchy", "l_v2x_polyarchy", 
                                                "s_lpec", "l_lpec", "yrs")) 

q10 <- q10 + scale_x_discrete(labels=c("Regime change", "Interstate war", "Intrastate war",
                                       "Polyarchy score (change)", "Polyarchy score (lag)",
                                       "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: War. Regime Change, and Changes in Ethnic Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_6.pdf", plot = q10, width = 10, height = 8)

##############################################################################################
######################################################################################
## Robustness Checks:
##############################################################################################


######################################################################################
# Appendix Figure A2-a: recent win
recent_win <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                  'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                  'f10s11_v2clsocgrp'),
                          ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                          modname = c('ethnicCL_outcome0', 'ethnicCL_outcome1', 'ethnicCL_outcome2', 'ethnicCL_outcome3', 'ethnicCL_outcome4',
                                      'ethnicCL_outcome5', 'ethnicCL_outcome10'),
                          n_sim = 1000,
                          var = "recent_win",
                          data = data_cty)
# add labels
p1 <- recent_win$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Recent Wins and Changes in Ethnic Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_2a.pdf", plot = p1, width = 6.5, height = 4)

## Appendix Fig A2-b: recent losses
recent_loss <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                 'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                 'f10s11_v2clsocgrp'),
                         ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                 "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                         modname = c('ethnicCL_outcome0', 'ethnicCL_outcome1', 'ethnicCL_outcome2', 'ethnicCL_outcome3', 'ethnicCL_outcome4',
                                     'ethnicCL_outcome5', 'ethnicCL_outcome10'),
                         n_sim = 1000,
                         var = "recent_loss",
                         data = data_cty)
# add labels
p2 <- recent_loss$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                    "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Recent Losses and Changes in Ethnic Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_2b.pdf", plot = p2, width = 6.5, height = 4)

# Appendix Fig2A-c: recent draws
recent_draw <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                  'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                  'f10s11_v2clsocgrp'),
                          ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                          modname = c('ethnicCL_outcome0', 'ethnicCL_outcome1', 'ethnicCL_outcome2', 'ethnicCL_outcome3', 'ethnicCL_outcome4',
                                      'ethnicCL_outcome5', 'ethnicCL_outcome10'),
                          n_sim = 1000,
                          var = "recent_draw",
                          data = data_cty)
# add labels
p3 <- recent_draw$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Recent Draws and Changes in Ethnic Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_2c.pdf", plot = p3, width = 6.5, height = 4)

# Appendix Fig 2A-d: ongoing wars
ongoing_war <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                  'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                  'f10s11_v2clsocgrp'),
                          ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                          modname = c('ethnicCL_outcome0', 'ethnicCL_outcome1', 'ethnicCL_outcome2', 'ethnicCL_outcome3', 'ethnicCL_outcome4',
                                      'ethnicCL_outcome5', 'ethnicCL_outcome10'),
                          n_sim = 1000,
                          var = "ongoing_war",
                          data = data_cty)
# add labels
p4 <- ongoing_war$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Ongoing Wars and Changes in Ethnic Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_2d.pdf", plot = p4, width = 6.5, height = 4)


## Appendix figure A4
q3 <- coef_rb_plot(ongoing_war$Results, vars = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                                 "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs")) 

q3 <- q3 + scale_x_discrete(labels=c("New war", "Ongoing war", "Recent win", "Recent loss",
                                     "Recent draw","Polyarchy score (change)", "Polyarchy score (lag)",
                                     "Energy (change)", "Energy(lag)", "Year")) + 
       ggtitle("Coefficient Plot: War Outcomes and Changes in Ethnic Civil-liberty Equality") + 
       scale_color_discrete(name = "", 
                            labels = c("Current", "1-year",  "2-year",
                                       "3-year", "4-year", "5-year",
                                       "10-year"))
ggsave("figures/Appendix_Fig_4.pdf", plot = q3, width = 11, height = 9)




###########################################################################################################
### Appendix Figure 3A: war outcome and changes in women's civil liberties

# Fig 3A-a: recent win
recent_win <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                 'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                 'f10s11_civilliberty'),
                         ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                 "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                         modname = c('womenCL_outcome0', 'womenCL_outcome1', 'womenCL_outcome2', 'womenCL_outcome3', 'womenCL_outcome4',
                                     'womenCL_outcome5', 'womenCL_outcome10'),
                         n_sim = 1000,
                         var = "recent_win",
                         data = data_cty)
# add labels
p1 <- recent_win$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                    "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Recent Wins and Changes in Gender Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_3a.pdf", plot = p1, width = 6.5, height = 4)



## Appendix Fig 3A-b: recent losses
recent_loss <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                  'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                  'f10s11_civilliberty'),
                          ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                          modname = c('womenCL_outcome0', 'womenCL_outcome1', 'womenCL_outcome2', 'womenCL_outcome3', 'womenCL_outcome4',
                                      'womenCL_outcome5', 'womenCL_outcome10'),
                          n_sim = 1000,
                          var = "recent_loss",
                          data = data_cty)
# add labels
p2 <- recent_loss$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Recent Losses and Changes in Gender Civil-liberty Equality",
            x = "Average marginal effects") + theme(plot.title = element_text(hjust = 0.6))
ggsave("figures/Appendix_Fig_3b.pdf", plot = p2, width = 6.5, height = 4)

# Appendix Fig 3A-c: recent draws
recent_draw <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                  'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                  'f10s11_civilliberty'),
                          ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                          modname = c('womenCL_outcome0', 'womenCL_outcome1', 'womenCL_outcome2', 'womenCL_outcome3', 'womenCL_outcome4',
                                      'womenCL_outcome5', 'womenCL_outcome10'),
                          n_sim = 1000,
                          var = "recent_draw",
                          data = data_cty)
# add labels
p3 <- recent_draw$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Recent Draws and Changes in Gender Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_3c.pdf", plot = p3, width = 6.5, height = 4)

# Appendix Fig 3A-d: ongoing wars
ongoing_war <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                  'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                  'f10s11_civilliberty'),
                          ivs = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs"),
                          modname = c('womenCL_outcome0', 'womenCL_outcome1', 'womenCL_outcome2', 'womenCL_outcome3', 'womenCL_outcome4',
                                      'womenCL_outcome5', 'womenCL_outcome10'),
                          n_sim = 1000,
                          var = "ongoing_war",
                          data = data_cty)
# add labels
p4 <- ongoing_war$plot + scale_y_discrete(labels = c("10-year", "5-year","4-year", 
                                                     "3-year", "2-year", "1-year", "Current")) + 
       labs(title = "Ongoing Wars and Changes in Gender Civil-liberty Equality",
            x = "Average marginal effects")
ggsave("figures/Appendix_Fig_3d.pdf", plot = p4, width = 6.5, height = 4)


## Appendix figure A5
q4 <- coef_rb_plot(ongoing_war$Results, vars = c('new_war' , 'ongoing_war',"recent_win","recent_loss","recent_draw",
                                                 "s_v2x_polyarchy", "l_v2x_polyarchy", "s_lpec", "l_lpec", "yrs")) 

q4 <- q4 + scale_x_discrete(labels=c("New war", "Ongoing war", "Recent win", "Recent loss",
                                     "Recent draw","Polyarchy score (change)", "Polyarchy score (lag)",
                                     "Energy (change)", "Energy(lag)", "Year")) + 
       ggtitle("Coefficient Plot: War Outcomes and Changes in Gender Civil-liberty Equality") + 
       scale_color_discrete(name = "", 
                            labels = c("Current", "1-year",  "2-year",
                                       "3-year", "4-year", "5-year",
                                       "10-year"))
ggsave("figures/Appendix_Fig_5.pdf", plot = q4, width = 11, height = 9)
#################################################################################


## Appendix Figure A8

### robustness check: December 13, 2019
data_cty_sub <- data_cty %>% 
            dplyr::filter(year > 1988)

## create two dummies for war and svac interaction
data_cty_sub <- data_cty_sub %>% 
          dplyr::mutate(interwar_sv = ifelse(inter_warDummy == 1 & sv == 1, 1, 0),
                        interwar_nosv = ifelse(inter_warDummy == 1 & sv == 0, 1, 0),
                        intrawar_sv = ifelse(intra_warDummy == 1 & sv == 1, 1, 0),
                        intrawar_nosv = ifelse(intra_warDummy == 1 & sv == 0, 1, 0))

interaction <- FE_rb_plot(dvs = c('s_civilliberty', 'fs2_civilliberty', 'f2s3_civilliberty',
                                  'f3s4_civilliberty', 'f4s5_civilliberty', 'f5s6_civilliberty',
                                  'f10s11_civilliberty'),
                          ivs = c('interwar_sv' , 'interwar_nosv','intrawar_sv','intrawar_nosv',
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", 
                                  "s_lpec", "l_lpec", "yrs"),
                          modname = c('womenCL0', 'womenCL1', 'womenCL2', 'womenCL3', 'womenCL4',
                                      'womenCL5', 'womenCL10'),
                          n_sim = 1000,
                          var = "interwar_sv",
                          data = data_cty_sub)
summary(interaction$Fits[[7]])
# The order of the labels will be the order of the input
coef_interaction <- coef_rb_plot(interaction$Results, vars = c(
              'interwar_sv' , 'interwar_nosv','intrawar_sv','intrawar_nosv',
              "s_v2x_polyarchy", "l_v2x_polyarchy", 
              "s_lpec", "l_lpec", "yrs")) 
#add lables 
coef_interaction <- coef_interaction + scale_x_discrete(labels=c("Interstate war w/ sexual violence",
                                                                 "Interstate war w/o sexual violence",
                                                                 "Intrastate war w/ sexual violence",
                                                                 "Intrastate war w/o sexual violence",
                                               "Polyarchy score (change)", "Polyarchy score (lag)",
                                               "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: Interactive effects of War and Changes in Gender Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_8a.pdf", plot = coef_interaction, width = 12, height = 6.5)


interaction_eth <- FE_rb_plot(dvs = c('s_v2clsocgrp', 'fs2_v2clsocgrp', 'f2s3_v2clsocgrp',
                                  'f3s4_v2clsocgrp', 'f4s5_v2clsocgrp', 'f5s6_v2clsocgrp',
                                  'f10s11_v2clsocgrp'),
                          ivs = c('interwar_sv' , 'interwar_nosv','intrawar_sv','intrawar_nosv',
                                  "s_v2x_polyarchy", "l_v2x_polyarchy", 
                                  "s_lpec", "l_lpec", "yrs"),
                          modname = c('ethnicityCL0', 'ethnicityCL1', 'ethnicityCL2', 'ethnicityCL3', 'ethnicityCL4',
                                      'ethnicityCL5', 'ethnicityCL10'),
                          n_sim = 1000,
                          var = "interwar_sv",
                          data = data_cty_sub)

# The order of the labels will be the order of the input
coef_interaction_eth <- coef_rb_plot(interaction_eth$Results, vars = c(
                    'interwar_sv' , 'interwar_nosv','intrawar_sv','intrawar_nosv',
                    "s_v2x_polyarchy", "l_v2x_polyarchy", 
                    "s_lpec", "l_lpec", "yrs")) 
#add lables 
coef_interaction_eth <- coef_interaction_eth + scale_x_discrete(labels=c("Interstate war w/ sexual violence",
                                                                 "Interstate war w/o sexual violence",
                                                                 "Intrastate war w/ sexual violence",
                                                                 "Intrastate war w/o sexual violence",
                                                                 "Polyarchy score (change)", "Polyarchy score (lag)",
                                                                 "Energy (change)", "Energy(lag)", "Year")) + 
  ggtitle("Coefficient Plot: Interactive effects of War and Changes in Ethnic Civil-liberty Equality") + 
  scale_color_discrete(name = "", 
                       labels = c("Current", "1-year",  "2-year",
                                  "3-year", "4-year", "5-year",
                                  "10-year"))
ggsave("figures/Appendix_Fig_8b.pdf", plot = coef_interaction_eth, width = 12, height = 6.5)


########################################################################################
## Figure 4: Group-level Analysis
########################################################################################
# load the group-level data
load("df_group.RData")

#### use lmer4
# fit a multilevel model
# group level results: Fig 4
mlm_pwrrank <- MLMs_fit(dvs = c('s_status_pwrrank', 'fs2_status_pwrrank',
                                'f2s3_status_pwrrank','f3s4_status_pwrrank', 
                                'f4s5_status_pwrrank', 'f5s6_status_pwrrank',
                                'f10s11_status_pwrrank'),
                        ivs = c('ucdp_terr', 'ucdp_gov','incidence_terr_flag', 
                                'incidence_gov_flag', 'democ', 'lcpop', 'year',
                                'elec_constituent_dum', 'elec_executive_dum'),
                        modname = c('pwrrankCL0', 'pwrrankCL1', 'pwrrankCL2', 
                                    'pwrrankCL3', 'pwrrankCL4',
                                    'pwrrankCL5', 'pwrrankCL10'),
                        data = df_group)

## Figure 4a 
mlm_pwrrank_ucdpTERR <- MLM_maginal_plot(Fits = mlm_pwrrank, n_sim = 1000,
                                         var = "ucdp_terr",
                                         data = df_group)
ucdp_terr <- mlm_pwrrank_ucdpTERR$plot +  
  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                              "2-year", "1-year", "Current")) + 
  labs(title = "Territorial War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Figure_4a.eps", plot = ucdp_terr, width = 6.5, height = 3.5, dpi = 400, units = "in")


# Fig 4b: UCDP governmental 
mlm_pwrrank_ucdpGOV <- MLM_maginal_plot(Fits = mlm_pwrrank, n_sim = 1000,
                                        var = "ucdp_gov",
                                        data = df_group)
ucdp_gov <- mlm_pwrrank_ucdpGOV$plot +  
  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                              "2-year", "1-year", "Current")) + 
  labs(title = "Governmental War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Figure_4b.eps", plot = ucdp_gov, width = 6.5, height = 3.5, dpi = 400, units = "in")


## Fig 4c: 
mlm_pwrrank_eprTERR <- MLM_maginal_plot(Fits = mlm_pwrrank, n_sim = 1000,
                                        var = "incidence_terr_flag",
                                        data = df_group)
epr_terr_flag <- mlm_pwrrank_eprTERR$plot +  
  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                              "2-year", "1-year", "Current")) + 
  labs(title = "Group's Territorial War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Figure_4c.eps", plot = epr_terr_flag, width = 6.5, height = 3.5, dpi = 400, units = "in")


## Fig 4d

mlm_pwrrank_eprGOV <- MLM_maginal_plot(Fits = mlm_pwrrank, n_sim = 1000,
                                       var = "incidence_gov_flag",
                                       data = df_group)
epr_gov_flag <- mlm_pwrrank_eprGOV$plot +  
  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                              "2-year", "1-year", "Current")) + 
  labs(title = "Group's Governmental War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Figure_4d.eps", plot = epr_gov_flag, width = 6.5, height = 3.5, dpi = 400, units = "in")


####################### Appendix Table A10
#To appendix table
screenreg(mlm_pwrrank, digits = 4)
texreg(mlm_pwrrank, file = "Appendix_TableA1.tex", digits = 4, label = "tab:ethnicpwrrank",
       custom.model.names = c("Current", "1-year", "2-year", "3-year", "4-year","5-year", "10-year"),
       custom.coef.names = c("Intercept", "Territorial war within country", "Governmental war within country",
                             "Territorial war with group", "Governmental war with group", "Democracy", "Country Population (log)",
                             "Year", "Constituent election", "Executive election"),
       caption.above = TRUE, caption = "War type and change in group's ethnic status ranking (mixed-effect model)")


############################################################################
####### Appendix Figure A7
############################################################################
# Appendix: group level
## using fixed effect OLS model

# Appendix Figure A7-a

FE_ethic <-  FEgroup_rb_plot(dvs = c('s_status_pwrrank', 'fs2_status_pwrrank',
                                     'f2s3_status_pwrrank','f3s4_status_pwrrank', 
                                     'f4s5_status_pwrrank', 'f5s6_status_pwrrank',
                                     'f10s11_status_pwrrank'),
                             ivs = c('ucdp_terr', 'ucdp_gov','incidence_terr_flag', 
                                     'incidence_gov_flag', 'democ', 'lcpop', 'yrs'),
                             modname = c('pwrrankCL0', 'pwrrankCL1', 'pwrrankCL2', 
                                         'pwrrankCL3', 'pwrrankCL4',
                                         'pwrrankCL5', 'pwrrankCL10'),
                             n_sim = 1000, 
                             var = "ucdp_terr",
                             data = df_group)
fe_ucdp_terr <- FE_ethic$plot +  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                                                             "2-year", "1-year", "Current")) + 
  labs(title = "Territorial War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Appendix_Fig_7a.pdf", plot = fe_ucdp_terr, width = 7, height = 6)

# Appendix Figure A7-b
FE_ethic <-  FEgroup_rb_plot(dvs = c('s_status_pwrrank', 'fs2_status_pwrrank',
                                     'f2s3_status_pwrrank','f3s4_status_pwrrank', 
                                     'f4s5_status_pwrrank', 'f5s6_status_pwrrank',
                                     'f10s11_status_pwrrank'),
                             ivs = c('ucdp_terr', 'ucdp_gov','incidence_terr_flag', 
                                     'incidence_gov_flag', 'democ', 'lcpop', 'yrs'),
                             modname = c('pwrrankCL0', 'pwrrankCL1', 'pwrrankCL2', 
                                         'pwrrankCL3', 'pwrrankCL4',
                                         'pwrrankCL5', 'pwrrankCL10'),
                             n_sim = 1000, 
                             var = "ucdp_gov",
                             data = df_group)
fe_ucdp_gov <- FE_ethic$plot +  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                                                            "2-year", "1-year", "Current")) + 
  labs(title = "Governmental War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Appendix_Fig_7b.pdf", plot = fe_ucdp_gov, width = 7, height = 6)


# Appendix Figure A7-c
FE_ethic <-  FEgroup_rb_plot(dvs = c('s_status_pwrrank', 'fs2_status_pwrrank',
                                     'f2s3_status_pwrrank','f3s4_status_pwrrank', 
                                     'f4s5_status_pwrrank', 'f5s6_status_pwrrank',
                                     'f10s11_status_pwrrank'),
                             ivs = c('ucdp_terr', 'ucdp_gov','incidence_terr_flag', 
                                     'incidence_gov_flag', 'democ', 'lcpop', 'yrs'),
                             modname = c('pwrrankCL0', 'pwrrankCL1', 'pwrrankCL2', 
                                         'pwrrankCL3', 'pwrrankCL4',
                                         'pwrrankCL5', 'pwrrankCL10'),
                             n_sim = 1000, 
                             var = "incidence_terr_flag",
                             data = df_group)
fe_epr_terr <- FE_ethic$plot +  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                                                            "2-year", "1-year", "Current")) + 
  labs(title = "Group's Territorial War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Appendix_Fig_7c.pdf", plot = fe_epr_terr, width = 7, height = 6)

# Appendix Figure A7-d
### use EPR 
FE_ethic <-  FEgroup_rb_plot(dvs = c('s_status_pwrrank', 'fs2_status_pwrrank',
                                     'f2s3_status_pwrrank','f3s4_status_pwrrank', 
                                     'f4s5_status_pwrrank', 'f5s6_status_pwrrank',
                                     'f10s11_status_pwrrank'),
                             ivs = c('ucdp_terr', 'ucdp_gov','incidence_terr_flag', 
                                     'incidence_gov_flag', 'democ', 'lcpop', 'yrs'),
                             modname = c('pwrrankCL0', 'pwrrankCL1', 'pwrrankCL2', 
                                         'pwrrankCL3', 'pwrrankCL4',
                                         'pwrrankCL5', 'pwrrankCL10'),
                             n_sim = 1000, 
                             var = "incidence_gov_flag",
                             data = df_group)
fe_epr_gov <- FE_ethic$plot +  scale_y_discrete(labels = c("10-year","5-year","4-year", "3-year",
                                                           "2-year", "1-year", "Current")) + 
  labs(title = "Group's Governmental War and Change in Group Status Ranking",
       x = "Average marginal effects")
ggsave("figures/Appendix_Fig_7d.pdf", plot = fe_epr_gov, width = 7, height = 6)

############# End of the replication code



